# -*- coding: utf-8 -*-
#libraries and modules


import math
import numpy as np
import scipy as sp
import scipy.stats.stats as st
import pandas as pd
import matplotlib.pyplot as plt
import time

import HBtoMCstrength as HBMC
import Limit_Equilibrium as LE
import SlipSurface as SS
import Slices as Sl



#math/trigo functions
tan=math.tan
sin=math.sin
atan=math.atan
rad=math.radians
pow=math.pow
exp=math.exp


def generate_HB_FE_list(hl_list,B_list,GSI_list,D,mi_list,sig_ci_list,uw):
    #DESCRIPTION:
    #generate Hoek-Brown failure envelopes for all sets of material properties

    #INPUT
    #hl_list;       numpy array of floats;  length of slope base, in m
    #B_list;        numpy array of floats;  slope gradient, in degrees
    #GSI_list;      numpy array of floats;  Geological Strength Index
    #D;             float; Disturbance factor
    #mi_list;       numpy array of floats; material constant for intact rock
    #sig_ci_list;   numpy array of floats; unconfined compressive strength (UCS) of intact rock, in MPa
    #uw;            float; unit weight of rock, in MN/m^3

    #OUTPUT
    #sig_t_list;    numpy array of floats; tensile strength of rock mass
    #sign_tau_list; list of UnivariateSpline objects, representing 
    

    sign_tau_list=[]
    sig_t_list=[]
    for im in range(len(GSI_list)):        
        maxH=round(hl_list.max()*tan(rad(B_list.max())),0)
        sign_tau, sig_t=HBMC.HoekBrownFE_spline(GSI_list[im],D,mi_list[im],sig_ci_list[im],uw,maxH)
        sign_tau_list.append(sign_tau)
        sig_t_list.append(sig_t)
    sig_t_list=np.asarray(sig_t_list)
    return sig_t_list, sign_tau_list


def slope_pts(H,B,xarr):
    #DESCRIPTION: solves for y coordinates of slope surface given the x-coordinates and slope height and gradient
    xarr=np.asarray(xarr)
    zero=np.zeros(len(xarr))
    crest=H*np.ones(len(xarr))
    y=np.where(xarr>0,xarr*tan(rad(B)),zero)
    y=np.where(xarr>H/tan(rad(B)),crest,y)
    return y

def print_Slide_info(H,B,num_slip,min_alpha,delta_B_max_alpha,minx1c,maxx2c):
    
    print "*********SLIDE INFO*********"
    print "\nBoundaries: "
    print [round(minx1c*H/sin(rad(B)),3),0], [0,0], [round(H/tan(rad(B)),3),round(H,3)],[round(H/tan(rad(B))+maxx2c*H/sin(rad(B)),3),round(H,3)] 
    print [round(H/tan(rad(B))+maxx2c*H/sin(rad(B)),3),round(-H/2.,3)],[round(minx1c*H/sin(rad(B)),3),round(-H/2.,3)]

    print "\nslipe surfaces"
    print " number of slip surfaces: ", num_slip
    print " slip surface limits:"
    print "     x1: ", [0, round(0.5*(H/tan(rad(B))),3)]
    print "     x2: ", [round(H/tan(rad(B)),3), round(H/tan(rad(B))+0.5*(H/sin(rad(B))),3)]
    print "     min exit angle: ", min_alpha
    print "     max exit angle: ", B-delta_B_max_alpha

    print "\nnumber of slices: ", num_slices

    print "\nmaterial properties: "
    print " GSI: ", GSI
    print " sigci: ", sig_ci
    print " mi: ", mi
    print " D: ", D
    print " UW: ", uw

    print "external forces"
    print " ru: ", ru
    print " kh: ", kh

    print "\ncheck m-alpha: ", malpha_check
    print "check negative normal force: ", neg_sign_check

    print "\n*********SLIDE INFO*********"
    

def plot_FS_results(FS,maxFS,h,k,R,x1,x2,H,B,minx1c,maxx2c,print_arc_centers):
    scaled_z = FS/maxFS
    colors = plt.cm.gist_rainbow(scaled_z)
    for g in np.argsort(-FS):
        xarc=np.linspace(x1[g],x2[g],50)
        yarc=k[g]-np.sqrt(R[g]**2-(xarc-h[g])**2)
        plt.plot(xarc,yarc,'-',c=colors[g])
        if print_arc_centers:
            plt.plot([h[g]],[k[g]],'+',c=colors[g])
    plt.plot(xarc,yarc,'k-')
    if print_arc_centers:
        plt.plot([h[g]],[k[g]],'ko')
    xarr=np.linspace(minx1c*H/tan(rad(B)),H/tan(rad(B))+maxx2c*H/sin(rad(B)),100)
    plt.plot(xarr, slope_pts(H,B,xarr),'k-',lw=1.5)
    
    return 

#############################
#                           #
#    USER INPUTS            #
#                           #
#############################


#random seed
np.random.seed(102345)
    
#1. SLOPE GEOMETRY
hl_list=np.asarray([30.,50.,70.,90.,110.])
B_list=np.asarray([15.,30.,45.,60.,75.])

#2. DYNAMIC FORCES A
kh=0.3              #seismic coefficient (horizontal)
ru=0.3              #pore pressure coefficient



#4. GENERATING SLICES
num_slices=10           #input number of slices for each slip surface
                
#5. SLIP SURFACE
fixed_surface=0         #1 for fixed, single value; 0 for multiple, random values
if fixed_surface:       #input single value for each parameter 
    x1=np.asarray([20.977])
    x2=np.asarray([59.296])
    h=np.asarray([-198.738])
    k=np.asarray([243.561])
    R=np.asarray([302.024])
    num_slip=1
else:
    num_slip=2000           #input number of slip surfaces to generate  
    min_alpha=-30           #minimum exit angle of slip surface 
    delta_B_max_alpha=10    #minimum difference between B maximum exit angle of slip surface
min_x1_c=-1             #factor of slope base length, <=0.5;                (+) within slope face; (0) toe; (-) slope base 
max_x1_c=0.5            #factor of slope base length, >=min_x1_c, <=0.5;    (+) within slope face; (0) toe; (-) slope base
min_x2_c=0              #factor of slope face length, >=0;                  (0) at slope crest; (+) beyond slope crest 
max_x2_c=0.5            #factor of slope face length, >=min_x2; 

#6. SLOPE STABILITY ANALYSIS CHECKS
malpha_check=0
neg_sign_check=1
converge_check=1

#7. FS OUTPUTS
full_output_FS=0     #1 full output of slice values will be returned, otherwise, only FS value is returned.
full_output_Slices=0
plot_deterministic=0
print_arc_centers=0
maxFS=6.0




#8. MATERIAL PROPERTIES
fixed_materials=0    #1 for fixed, single value; 0 for multiple, random values
if fixed_materials:  #input single value for each property 
    lith_list=['fixed'] 

    #geologic strength index
    GSI_list=np.asarray([49.])

    #material constant
    mi_list=np.asarray([17.])

    #unconfined compressive strength, MPa
    sig_ci_list=np.asarray([15.7])

    #unit weight, MN/m3
    uw=25*10**-3

    #disturbance factor
    D=0.7

else:
    lith_list=['CC','FC','LS','UM']     #input list of lithologic groups in external csv files
    UCS_rf_list=[0.5, 0.3, 0.65, 0.75]
    num_mats=500                        #input number of sets of material properties for probabilistic analysis

    triang_headers=['Li','lft','mde','rht'] #input headers for *_triang.csv files

    #reading properties from external csv files
    DF_GSI=pd.read_csv('GSI_triang.csv',sep=',',names=triang_headers,index_col=[0])
    DF_UCS=pd.read_csv('UCS_triang.csv',sep=',',names=triang_headers,index_col=[0])
    DF_mi=pd.read_csv('mi_triang.csv',sep=',',names=triang_headers,index_col=[0])
    DF_uw=pd.read_csv('uw_constant.csv',sep=',',names=['Li','kval'],index_col=[0])

    
    #generating random variables from triangular distribution
    list_GSI_list=[]
    list_UCS_list=[]
    list_mi_list=[]
    list_uw=[]
    for iLith in range(len(lith_list)):
        
        GSI=DF_GSI[(DF_GSI.index==lith_list[iLith])]
        UCS=DF_UCS[(DF_UCS.index==lith_list[iLith])]
        mi=DF_mi[(DF_mi.index==lith_list[iLith])]
        uw=DF_uw[(DF_uw.index==lith_list[iLith])]

        list_GSI_list.append(np.round(np.random.triangular(GSI.lft,GSI.mde,GSI.rht,num_mats),1))
        list_UCS_list.append(np.round(np.random.triangular(UCS.lft,UCS.mde,UCS.rht,num_mats),1))
        list_mi_list.append(np.round(np.random.triangular(mi.lft,mi.mde,mi.rht,num_mats),1))
        list_uw.append(np.round(uw.kval.values/1000.,3)[0])
   
    #disturbance factor
    D=0.7



#############################
#                           #
#    MAIN                   #
#                           #
#############################

##fig,ax=plt.subplots(nrows=2,ncols=4,sharex=True, sharey=True)


t0=time.time()
#0. Looping through each lithologic group
for iLith in range(len(lith_list)):
    if iLith!=3:continue
    if lith_list!=['fixed']:
        GSI_list=list_GSI_list[iLith]
        sig_ci_list=list_UCS_list[iLith]*UCS_rf_list[iLith]
        mi_list=list_mi_list[iLith]
        uw=round(list_uw[iLith],3)
    
##    curax=ax[0,iLith]
##    curax.hist(mi_list)
##
##    curax=ax[1,iLith]
##    curax.hist(mi_list*0.5)
##
##    continue

    


    print "\n",lith_list[iLith]
    #print GSI_list
    #print sig_ci_list
    #print mi_list
    #print uw

    #continue
    

    #1. Generate list of failure envelopes (sig_t_list, sign_tau_list) corresponding to the list of material properties
    sig_t_list, sign_tau_list=generate_HB_FE_list(hl_list,B_list,GSI_list,D,mi_list,sig_ci_list,uw)

##    sig_t_list2, sign_tau_list2=generate_HB_FE_list(hl_list,B_list,GSI_list,D,mi_list,sig_ci_list*0.5,uw)
##
##
##    for im in range(len(sign_tau_list)):
##        #setting current failure envelope
##        sign_tau=sign_tau_list[im]
##        sign_tau2=sign_tau_list2[im]
##
##        x=np.linspace(0,3,1000)
##        plt.plot(x,sign_tau(x),'b-')
##        plt.plot(x,sign_tau2(x),'r-')
##        plt.xlim(0,0.1)
##        plt.show()
##    continue


    #2. Initialize figure
    fig=plt.figure(facecolor=None, edgecolor=None)
    plt.title(lith_list[iLith])
    plt.axis('off')

    #3. looping through each hl
    for ihl in range(len(hl_list)):

        #4. looping through each B
        for iB in range(len(B_list)):
            

            #5. defining current H,B and corresponding figure subplot
            B=B_list[len(B_list)-iB-1]
            H=round(hl_list[ihl]*tan(rad(B)),2)
            print "\n",hl_list[ihl],B
            #creating subplot
            if B==B_list[-1]:
                curax=fig.add_subplot(len(B_list), len(hl_list), (len(hl_list)*iB+1)+ihl)
                ax_sharex=curax
            else:
                curax=fig.add_subplot(len(B_list), len(hl_list), (len(hl_list)*iB+1)+ihl, sharex=ax_sharex)
            #defining x-axis limit    
            xax_min,xax_max=0-hl_list[ihl],hl_list[ihl] + hl_list[ihl]*tan(rad(B_list.max()))/sin(rad(B_list.max())) 
            #setting and formatting current axis
            plt.sca(curax)
            if B==B_list.min():
                curax.set_xlabel(str(hl_list[ihl]))
            else:
                curax.set_xticklabels([])
            if ihl==0:
                curax.set_ylabel(str(B))
            

            #6. GENERATING SLIP SURFACES FOR CURRENT H,B
            

            if fixed_surface==0:
                
               
                x1,x2,h,k,R,errors=SS.slip_surface_generation(H,B,num_slip, min_x1_c,max_x1_c,min_x2_c,max_x2_c, min_alpha, delta_B_max_alpha,num_slices)
                                                                        #print "Err105 (>2 slip surface - slope intsxns): ",len(np.where(errors==-105)[0])
                                                                        #print "Err114 (reverse curvature): ",len(np.where(errors==-114)[0])

                print "valid slices: ", len(np.where(errors==1)[0])
            
            
        

            
            #7. GENERATING SLICE VALUES FOR ALL SLIP SURFACES
            slice_values=[]
            for c in range(len(x1)):
                weight,alpha,bl,pwp,SF,seisF,tcF, slices=Sl.Slice_params(H,B,num_slices,h[c],k[c],R[c],x1[c],x2[c],uw,ru,kh, full_output_Slices)
                slice_values.append([weight,alpha,bl,pwp,SF,seisF,tcF, slices])
            
            #initializong list of FS results for current H,B
            HB_FS=[]
            HB_ch=[]
            HB_ck=[]
            HB_cR=[]
            HB_cx1=[]
            HB_cx2=[]
            #8. looping through each failure envelope 
            for im in range(len(sign_tau_list)):
                #setting current failure envelope
                sign_tau=sign_tau_list[im]
                sig_t=sig_t_list[im]
                if im%10==0: 
                    print im,                                               

                #initializong list of FS results for current H,B and material property
                FS_bis_list=[]
                ch=[]
                ck=[]
                cR=[]
                cx1=[]
                cx2=[]

                #9. looping through each set of slice values (of the corresponding slip surfaces)
                errorsFS=[]
                for c in range(len(slice_values)):
                    #10. GETTING CURRENT SLICE VALUES
                    s_weight=slice_values[c][0]
                    s_alpha=slice_values[c][1]
                    s_bl=slice_values[c][2]
                    s_pwp=slice_values[c][3]
                    s_SF=slice_values[c][4]
                    s_seisF=slice_values[c][5]
                    s_tcF=slice_values[c][6]
                    s_slices=slice_values[c][7]
                                                                            #print s_slices

                    #11. SOLVING FS using Fellenius method
                    FS_fel, phi_fel, coh_fel, s_fel=LE.OrdinaryMethod(s_weight,s_alpha,s_pwp,s_bl,s_SF,s_seisF,s_tcF,sign_tau, sig_t,neg_sign_check,full_output_FS)
                                                                            #print s_fel
                                                                            #print "FS_fel: ",FS_fel
                    if FS_fel<0:
                        errorsFS.append(FS_fel)
                        continue

                    #12. SOLVING FS using Simplified Bishop method
                    FS_bis, s_bis=LE.SimplifiedBishop(FS_fel,phi_fel,coh_fel,s_weight,s_bl,s_alpha,s_pwp,s_SF,s_tcF,s_seisF,sign_tau, sig_t,malpha_check,neg_sign_check,converge_check,full_output_FS)
                    if FS_bis>0:
                        FS_bis_list.append(FS_bis)
                        ch.append(h[c])
                        ck.append(k[c])
                        cR.append(R[c])
                        cx1.append(x1[c])
                        cx2.append(x2[c])
                        errorsFS.append(1)
                    else:
                        errorsFS.append(FS_bis)
                                                                            #print s_bis
                errorsFS=np.asarray(errorsFS)

                                                                            #print "Err111 (solution did not converge): ",len(np.where(errorsFS==-111)[0])
                                                                            #print "Err112 (malpha<0.2): ",len(np.where(errorsFS==-112)[0])
                                                                            #print "Err120 (neg eff normal stress): ",len(np.where(errorsFS==-120)[0])
                                                                            #print "valid FS: ", len(np.where(errorsFS==1)[0])

                
                #13. Deterministic results for current H,B and material property (Global Minimum)
                FS_bis_list=np.asarray(FS_bis_list)
                ch=np.asarray(ch)
                ck=np.asarray(ck)
                cR=np.asarray(cR)
                cx1=np.asarray(cx1)
                cx2=np.asarray(cx2)
                exit_pt=cx1*100/(H/tan(rad(B)))
                if plot_deterministic:
                    plt.figure(100)
                    plot_FS_results(FS_bis_list,maxFS,ch,ck,cR,cx1,cx2,H,B,min_x1_c,max_x2_c,print_arc_centers)
                    print ""
                    print "min FS: ", round(FS_bis.min(),2)
                    print "mean FS: ", round(FS_bis.mean(),2)
                    print "max FS: ", round(FS_bis.max(),2)
                    print "% FS<1: ",round(st.percentileofscore(FS_bis,1.0),2)
                    print "toe exits: ", np.mean(exit_pt)
                    plt.show()
                HB_FS.append(FS_bis_list.min())
                HB_ch.append(ch[np.argmin(FS_bis_list)])
                HB_ck.append(ck[np.argmin(FS_bis_list)])
                HB_cR.append(cR[np.argmin(FS_bis_list)])
                HB_cx1.append(cx1[np.argmin(FS_bis_list)])
                HB_cx2.append(cx2[np.argmin(FS_bis_list)])
                

            #14. Probabilistic results for current H,B and all material properties (Overall slope)
            HB_FS=np.asarray(HB_FS)
            HB_ch=np.asarray(HB_ch)
            HB_ck=np.asarray(HB_ck)
            HB_cR=np.asarray(HB_cR)
            HB_cx1=np.asarray(HB_cx1)
            HB_cx2=np.asarray(HB_cx2)
            HB_exit_pt=HB_cx1*100/(H/tan(rad(B)))
            plot_FS_results(HB_FS,maxFS,HB_ch,HB_ck,HB_cR,HB_cx1,HB_cx2,H,B,min_x1_c,max_x2_c,print_arc_centers)
            print ""
##            print "min FS: ", round(HB_FS.min(),2)
##            print "mean FS: ", round(HB_FS.mean(),2)
##            print "max FS: ", round(HB_FS.max(),2)
##            print "% FS<1: ",round(st.percentileofscore(HB_FS,1.0),2)
##            print "toe exits: ", np.mean(HB_exit_pt)

            print round(HB_FS.min(),2), round(HB_FS.mean(),2), round(HB_FS.max(),2), round(st.percentileofscore(HB_FS,1.0),2), np.mean(HB_exit_pt)

            
            plt.xlim(xax_min, xax_max)
            plt.axis('equal')
            print "time: ",round(time.time()-t0,2)
            
    #fig.tight_layout()
plt.show()

        


       
        
#references:
#Duncan JM and Wright SG (2005). Soil Strength and Slope Stability. John Wiley & Sons, Inc., Hoboken, New Jersey
#Hoek, E., Carranza-torres, C., & Corkum, B. (2002). Hoek-brown failure criterion – 2002 edition. NARMS-TAC (pp. 267-273).
#Rocscience. Slide online documentation. http://www.rocscience.com/help/slide/webhelp/Slide.htm

        

            
            
            

                

    











    
